PPFD is photosynthetic photon flux density in mol/m2/s Irradiance is the radiant flux received by a surface per unit of area in W/m2 Energy content in the PAR waveband = 2.35*10^5 J/mol
# irr = irradiance
PPFD<-function(irr){ # here the function called PPFD used the irradiance "irr" as an input
result<-irr*1/(2.35*10^5) # result represents the output (that is the PPFD)
return(result) # you need to use return to be sure the result will be printed in your R console
}
PPFD(500) # as PPFD function was created now you can put any input (irradiance)
## [1] 0.00212766
# if irr = 500 W/m2, PPFD is equal to 2100 micromol/m2/s
values<-seq(0,1000,by=100)
# values is a vector including value of irradiance fromm 0 to 1000 with a increment of 100
ggplot(data=NULL,aes(x=values,y=PPFD(values)))+geom_point() +
ylab("PPFD (mol/m2/s)") +
xlab("Irradiance (W/m2)") +
ggtitle("PPFD as a function of irradiance")
Planck’s law describes the spectral density of electromagnetic radiation emitted by a black body in thermal equilibrium at a given temperature T
h<-6.626*10^-34 # Planck's constant (J*s)
c<-2.997925*10^8 # speed of light (m/s) or 10^17 nm/s
k<-1.381*10^-23 # Boltzmann's constant (J/K)
planck<-function(wl,temp){ # wl = wavelength (m), T = temperature (K)
num<-2*pi*h*c^2
den<-wl^5*exp((h*c/k+wl+temp)-1)
result<-num/den
return(result)
}
h=6.62606896e-34
c=2.99792458e8
k=1.3806504e-23
lambda = seq(1e-8,1e-4, by=1e-8)
T = 6000
y=(2*pi*h*(c^2))/(lambda^5)*(1./(exp((h*c)/(lambda*k*T))-1))
T2 = 288
y2=(2*pi*h*(c^2))/(lambda^5)*(1./(exp((h*c)/(lambda*k*T2))-1))
plot(1e6*lambda,1e-6*1e-6*y, xlim=c(0.0,5.0),ylim=c(0.0,100.0), xlab='Wavelength (um)', ylab='spectral emittance (MW/m2/um)')
# par(new=T)
# plot(1e6*lambda,1e-6*1e-6*y2, xlim=c(0.0,10.0), ylim=c(0.0,1/10000.0), xlab='', ylab='', axes=F)
# par(new=F)
plot(1e6*lambda,1e-6*1e-6*y2, xlim=c(0.0,50.0), ylim=c(0.0,3*1e-5), xlab='Wavelength (um)', ylab='spectral emittance (MW/m2/um)')
lambda<-0.505 # micro m (or 505 nm)
# area pupil = 2.83*10^-5 m2
# 4.00*10^-11 W/m2 is the intensity for dark- adapted vision for human
# intensity per area = 4.00*10^-11 * 2.83*10^-5 = 1.13 *10^-15 J/s (reminder 1 watt = 1 J/s)
intensity<-1.13 *10^-15
# energy photon = 3.94*10^-19 J
# E = hc/lambda (eV)
hc<-1.24 *10^-6 # eV/ m
hc<-1.24 # eV/micro m
E<-hc/lambda
# one eV = 1.602 *10^-19 J
result<-1.602 *10^-19*E # result in J
number_photons<-intensity/result
numberphotons<-function(lambda) { # see PPFD function comments to understand how work the function. Input = lambda
intensity<-1.13 *10^-15
hc<-1.24
E<-hc/lambda
denom<-1.602 *10^-19*E
result<-intensity/denom
return(result)
}
lambda_values<-seq(400,800,by=50)
ggplot(data=NULL,aes(x=lambda_values,y=numberphotons(lambda_values))) + ### ggplot use a dataset (here null because I used a vector of values, aes is to define x and y)
geom_point() + # I used geom_point because I want points in the graph, for example if you want a line you should use geom_line
ylab("number of photons") + # specify y axis legend
xlab("wavelength (nm)") + # specify x axis legend
ggtitle("Dark vision: Number of photons as a function of wavelenght") # specify the title
### ggplot use a dataset (here null because I used a vector of values, aes is to define x and y)
#pick one CO2 value (this is from Zhu's paper Fig.3)
#present 400ppm
CO2<-0.0004
#past #220ppm,parts-per-million
CO2<-0.00022
#future #700ppm
CO2<-0.0007
#this is O2 consentration
O2<-0.21
cat("CO2 concentration:",CO2,"%","\n")
## CO2 concentration: 7e-04 %
cat("O2 concentration:",O2,"%","\n")
## O2 concentration: 0.21 %
# step 1:get into plant (physically diffusion)
ratio1<-O2/CO2
cat("The O2/CO2 ratio of physically diffusion:",ratio1,"%","\n")
## The O2/CO2 ratio of physically diffusion: 300 %
# step 2:water dissolving
alpha<-25 #This is the CO2/O2 Solubility ratio at 25 degree
ratio2<-ratio1*1/alpha
cat("The O2/CO2 ratio of water dissolving:",ratio2,"%","\n")
## The O2/CO2 ratio of water dissolving: 12 %
# step 3:
tau<-75 #This is the CO2/O2 Specificity ratio
ratio3<-tau/ratio2
cat("The O2/CO2 ratio of water dissolving:",ratio3,"%","\n")
## The O2/CO2 ratio of water dissolving: 6.25 %
#results:
#past ratio3:1.964286
#future ratio3:6.25
# equations are from paper: Improved temperature response functions for models of Rubisco-limited photosynthesis
###this is the function of tau with temperature (equation 6 in the paper)
#This function will call another function:parameter(Tk)
#inputs: c:scaling constant (dimensionless);
# Ha:energy of activation (kJ mol-1);
R=8.31445 #molar gas constant, M2kgs-2K-1mol-1
Temp.f = function(c,Ha){
parameter = function(Tk){
#Tk is leaf temeprature
exp(c-Ha*10^3/(R*(Tk+273.15)))
}
parameter
}
#This is the function of leaf temeprature
#all c and Ha value are from table 1 in the paper
#these four variables are functions with input Tk
Vc.max = Temp.f(26.35,65.33) #maximum RuBP saturated rate of carboxylation (mmol m-2 s-1);
Vo.max = Temp.f(22.98,60.11) #maximum RuBP saturated rate of oxygenation (mmol m-2 s-1);
Ko = Temp.f(20.3,36.38) # Michaelis constant for O2 (mmol mol-1);
Kc = Temp.f(38.05,79.43 )# Michaelis constant for CO2 (mmol mol-1);
#efficiency
#this is from Zhu's paper
max.conversion=0.487 #we get this value from page 154
#input: CO2 concentration
# O2 concentration
# tempreature
#output:Decrease value by photorespiration
# Efficiency
efficiency = function (CO2,O2,temp){
tau = Vc.max(temp)*Ko(temp)*1000/(Kc(temp)*Vo.max(temp))
phi = O2/(CO2*tau) #this is the ratio of RuBP oxygenation to carboxylation
dpr = 1-3*(1-0.5*phi)/(3+3.5*phi) #decrease in caused by photorespiration
conversion=max.conversion-dpr#energy conversion efficiency
list(dpr,conversion)
}
#prepare for plot
#temperature range
leaf.temp=c(5:40)
#CO2 concentration
CO2=0.00038
#O2 concentration
O2=0.21
#efficiency
effici = efficiency(CO2,O2,leaf.temp)
dpr=effici[[1]]
conversion=effici[[2]]
#plot-------------------------------------
#plot four intermediate variables
plot(leaf.temp,Vc.max(leaf.temp),type="l",xlab="Temperature",ylab="Maximum RuBP saturated rate of carboxylation (mmol m-2 s-1)")
plot(leaf.temp,Vo.max(leaf.temp),type="l",xlab="Temperature",ylab="Maximum RuBP saturated rate of oxygenation (mmol m-2 s-1)")
plot(leaf.temp,Ko(leaf.temp),type="l",xlab="Temperature",ylab="Michaelis constant for O2 (mmol mol-1)")
plot(leaf.temp,Kc(leaf.temp),type="l",xlab="Temperature",ylab="Michaelis constant for CO2 (mmol mol-1)")
#Decrease value by photorespiration channges with temperature
#dimensionless
plot(leaf.temp,dpr,type="l",xlab="Temperature",ylab="Decrease by Photorespiration")
#energy conversion efficiency channges with temperature
#dimensionless
plot(leaf.temp,conversion,type="l",xlab="Temperature",ylab="Energy Conversion Efficiency")
library(readr)
#tab<-read_csv("https://raw.githubusercontent.com/Agron508/Homework/master/SS-110_1010_12918noon.csv")
tab<-read.csv("https://raw.githubusercontent.com/Agron508/Homework/master/SS-110_1010_12918noon.csv")
library(ggplot2)
Some cleaning
data<-tab[2:482,] # select rows
colnames(data)[1]<-"waveL"
data$waveL<-as.numeric(as.vector(data$waveL)) # to get a x_continous scale, data need to be numeric not factor
Plot the Irradiance as a function of wavelenght (nm)
time<-data[,84]# here I just pick a column randomly
class(time)
## [1] "factor"
# time has to be converted into numeric
as.numeric.factor <- function(x) {as.numeric(levels(x))[x]} # function to convert all the level of the factor into numeric without losing information
time<-as.numeric.factor(time)
## Warning in as.numeric.factor(time): NAs introduced by coercion
ggplot(data=data,aes(x=waveL,y=time))+ # see part 4.3 to find some comments about ggplot
geom_line() +
scale_x_continuous(breaks = round(seq(min(data$waveL), max(data$waveL), by = 100),1)) +
scale_y_continuous(breaks = round(seq(min(time), max(time), by = 0.005),1)) +
xlab("wavelength (nm)") +
ylab(expression(paste("Spectral Irradiance (W/m2/",mu,"m)")))
St = total incident solar radiation during growing season of which 48.7% is photosynthetically active (MJ/m2) Ei = light interception efficiency (between 50-75%) Ec = energy conversion efficiency Ep = partitioning efficiency or harvest index (here ignored) Yp = yield potential (MJ/m2)
Seed_ec = seed energy content (MJ/kg) = 22.2 in average Leaf_ec = leaf energy content (MJ/kg) = 19 in average Stem_ec = stem energy content (MJ/kg) = 17 in average Total_ec = average energy content for the plant = (22.2+17+19)/3=19.4 (MJ/kg)
St=2000
Ei=0.6
Ec = 0.06 # at 35°C
Yp = St*0.478*Ei*Ec
Total_ec=19.4
biomass=Yp/Total_ec # biomass (carbon) of the whole plant
# the biomass is equal to 1.77 kg/m2
x<-c(0.5,0.6,0.7,0.75) # x = Ei
y<-c(0.05,0.1,0.2,0.3) # y = Ec
St<-2000
data<-expand.grid(x,y) # gives all the combinations of x and y
colnames(data)<-c("x","y")
data$z<-St*0.478*data$x*data$y # biomass in MJ/m2
library("plot3D")
scatter3D(data$x, data$y, data$z, phi = 0, bty = "g",
pch = 20, cex = 1.8,
ticktype = "detailed",
xlab = "light interception efficiency",
ylab ="energy conversion efficiency",
zlab = "yield potential (MJ/m2)",
main="Yield potential (MJ/m2), Ei and Ec")
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
x<-c(0.5,0.6,0.7,0.75) # x = Ei
y<-c(0.05,0.1,0.2,0.3) # y = Ec
St<-2000
data<-expand.grid(x,y) # gives all the combinations of x and y
colnames(data)<-c("x","y")
data$z<-St*0.478*data$x*data$y # biomass in MJ/m2
p <- plot_ly(data, x = ~x, y = ~y, z = ~z, color = ~z) %>%
add_markers() %>%
layout(scene = list(xaxis = list(title = 'Ei'),
yaxis = list(title = 'Ec'),
zaxis = list(title = 'Yp')))
p
## need to find a way to calculate C
cass <- function(result){
I <- 0.035 ## with the assump. that r = 3 and O = .21 equation (3) Irradiance
vc <- 98 ## constant taken form Fig 2 explaination it is a constant
C <- 270 ## micromol per mol
Rd <- result
## A = (1-I/C)vc-rd
A <- (1 -(I/C))*vc -Rd
return(A)
}